home *** CD-ROM | disk | FTP | other *** search
/ AGA Toolkit '97 / The AGA Toolkit '97.iso / miscellaneous / science / maths / calc / source / qmath.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-07  |  21.0 KB  |  1,256 lines
  1. /*
  2.  * Copyright (c) 1994 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Extended precision rational arithmetic primitive routines
  7.  */
  8.  
  9. #include "qmath.h"
  10.  
  11.  
  12. NUMBER _qzero_ =    { { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
  13. NUMBER _qone_ =        { { _oneval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
  14. static NUMBER _qtwo_ =    { { _twoval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
  15. static NUMBER _qten_ =    { { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
  16. NUMBER _qnegone_ =    { { _oneval_, 1, 1 }, { _oneval_, 1, 0 }, 1 };
  17. NUMBER _qonehalf_ =    { { _oneval_, 1, 0 }, { _twoval_, 1, 0 }, 1 };
  18.  
  19.  
  20. /*
  21.  * Create another copy of a number.
  22.  *    q2 = qcopy(q1);
  23.  */
  24. NUMBER *
  25. qcopy(q)
  26.     register NUMBER *q;
  27. {
  28.     register NUMBER *r;
  29.  
  30.     r = qalloc();
  31.     r->num.sign = q->num.sign;
  32.     if (!zisunit(q->num)) {
  33.         r->num.len = q->num.len;
  34.         r->num.v = alloc(r->num.len);
  35.         zcopyval(q->num, r->num);
  36.     }
  37.     if (!zisunit(q->den)) {
  38.         r->den.len = q->den.len;
  39.         r->den.v = alloc(r->den.len);
  40.         zcopyval(q->den, r->den);
  41.     }
  42.     return r;
  43. }
  44.  
  45.  
  46. /*
  47.  * Convert a number to a normal integer.
  48.  *    i = qtoi(q);
  49.  */
  50. long
  51. qtoi(q)
  52.     register NUMBER *q;
  53. {
  54.     long i;
  55.     ZVALUE res;
  56.  
  57.     if (qisint(q))
  58.         return ztoi(q->num);
  59.     zquo(q->num, q->den, &res);
  60.     i = ztoi(res);
  61.     zfree(res);
  62.     return i;
  63. }
  64.  
  65.  
  66. /*
  67.  * Convert a normal integer into a number.
  68.  *    q = itoq(i);
  69.  */
  70. NUMBER *
  71. itoq(i)
  72.     long i;
  73. {
  74.     register NUMBER *q;
  75.  
  76.     if ((i >= -1) && (i <= 10)) {
  77.         switch ((int) i) {
  78.             case 0: q = &_qzero_; break;
  79.             case 1: q = &_qone_; break;
  80.             case 2: q = &_qtwo_; break;
  81.             case 10: q = &_qten_; break;
  82.             case -1: q = &_qnegone_; break;
  83.             default: q = NULL;
  84.         }
  85.         if (q)
  86.             return qlink(q);
  87.     }
  88.     q = qalloc();
  89.     itoz(i, &q->num);
  90.     return q;
  91. }
  92.  
  93.  
  94. /*
  95.  * Convert a number to a normal unsigned integer.
  96.  *    u = qtou(q);
  97.  */
  98. FULL
  99. qtou(q)
  100.     register NUMBER *q;
  101. {
  102.     FULL i;
  103.     ZVALUE res;
  104.  
  105.     if (qisint(q))
  106.         return ztou(q->num);
  107.     zquo(q->num, q->den, &res);
  108.     i = ztou(res);
  109.     zfree(res);
  110.     return i;
  111. }
  112.  
  113.  
  114. /*
  115.  * Convert a normal unsigned integer into a number.
  116.  *    q = utoq(i);
  117.  */
  118. NUMBER *
  119. utoq(i)
  120.     FULL i;
  121. {
  122.     register NUMBER *q;
  123.  
  124.     if (i <= 10) {
  125.         switch ((int) i) {
  126.             case 0: q = &_qzero_; break;
  127.             case 1: q = &_qone_; break;
  128.             case 2: q = &_qtwo_; break;
  129.             case 10: q = &_qten_; break;
  130.             default: q = NULL;
  131.         }
  132.         if (q)
  133.             return qlink(q);
  134.     }
  135.     q = qalloc();
  136.     utoz(i, &q->num);
  137.     return q;
  138. }
  139.  
  140.  
  141. /*
  142.  * Create a number from the given integral numerator and denominator.
  143.  *    q = iitoq(inum, iden);
  144.  */
  145. NUMBER *
  146. iitoq(inum, iden)
  147.     long inum, iden;
  148. {
  149.     register NUMBER *q;
  150.     long d;
  151.     BOOL sign;
  152.  
  153.     if (iden == 0)
  154.         math_error("Division by zero");
  155.     if (inum == 0)
  156.         return qlink(&_qzero_);
  157.     sign = 0;
  158.     if (inum < 0) {
  159.         sign = 1;
  160.         inum = -inum;
  161.     }
  162.     if (iden < 0) {
  163.         sign = 1 - sign;
  164.         iden = -iden;
  165.     }
  166.     d = iigcd(inum, iden);
  167.     inum /= d;
  168.     iden /= d;
  169.     if (iden == 1)
  170.         return itoq(sign ? -inum : inum);
  171.     q = qalloc();
  172.     if (inum != 1)
  173.         itoz(inum, &q->num);
  174.     itoz(iden, &q->den);
  175.     q->num.sign = sign;
  176.     return q;
  177. }
  178.  
  179.  
  180. /*
  181.  * Add two numbers to each other.
  182.  *    q3 = qadd(q1, q2);
  183.  */
  184. NUMBER *
  185. qadd(q1, q2)
  186.     register NUMBER *q1, *q2;
  187. {
  188.     NUMBER *r;
  189.     ZVALUE t1, t2, temp, d1, d2, vpd1, upd1;
  190.  
  191.     if (qiszero(q1))
  192.         return qlink(q2);
  193.     if (qiszero(q2))
  194.         return qlink(q1);
  195.     r = qalloc();
  196.     /*
  197.      * If either number is an integer, then the result is easy.
  198.      */
  199.     if (qisint(q1) && qisint(q2)) {
  200.         zadd(q1->num, q2->num, &r->num);
  201.         return r;
  202.     }
  203.     if (qisint(q2)) {
  204.         zmul(q1->den, q2->num, &temp);
  205.         zadd(q1->num, temp, &r->num);
  206.         zfree(temp);
  207.         zcopy(q1->den, &r->den);
  208.         return r;
  209.     }
  210.     if (qisint(q1)) {
  211.         zmul(q2->den, q1->num, &temp);
  212.         zadd(q2->num, temp, &r->num);
  213.         zfree(temp);
  214.         zcopy(q2->den, &r->den);
  215.         return r;
  216.     }
  217.     /*
  218.      * Both arguments are true fractions, so we need more work.
  219.      * If the denominators are relatively prime, then the answer is the
  220.      * straightforward cross product result with no need for reduction.
  221.      */
  222.     zgcd(q1->den, q2->den, &d1);
  223.     if (zisunit(d1)) {
  224.         zfree(d1);
  225.         zmul(q1->num, q2->den, &t1);
  226.         zmul(q1->den, q2->num, &t2);
  227.         zadd(t1, t2, &r->num);
  228.         zfree(t1);
  229.         zfree(t2);
  230.         zmul(q1->den, q2->den, &r->den);
  231.         return r;
  232.     }
  233.     /*
  234.      * The calculation is now more complicated.
  235.      * See Knuth Vol 2 for details.
  236.      */
  237.     zquo(q2->den, d1, &vpd1);
  238.     zquo(q1->den, d1, &upd1);
  239.     zmul(q1->num, vpd1, &t1);
  240.     zmul(q2->num, upd1, &t2);
  241.     zadd(t1, t2, &temp);
  242.     zfree(t1);
  243.     zfree(t2);
  244.     zfree(vpd1);
  245.     zgcd(temp, d1, &d2);
  246.     zfree(d1);
  247.     if (zisunit(d2)) {
  248.         zfree(d2);
  249.         r->num = temp;
  250.         zmul(upd1, q2->den, &r->den);
  251.         zfree(upd1);
  252.         return r;
  253.     }
  254.     zquo(temp, d2, &r->num);
  255.     zfree(temp);
  256.     zquo(q2->den, d2, &temp);
  257.     zfree(d2);
  258.     zmul(temp, upd1, &r->den);
  259.     zfree(temp);
  260.     zfree(upd1);
  261.     return r;
  262. }
  263.  
  264.  
  265. /*
  266.  * Subtract one number from another.
  267.  *    q3 = qsub(q1, q2);
  268.  */
  269. NUMBER *
  270. qsub(q1, q2)
  271.     register NUMBER *q1, *q2;
  272. {
  273.     NUMBER *r;
  274.  
  275.     if (q1 == q2)
  276.         return qlink(&_qzero_);
  277.     if (qiszero(q2))
  278.         return qlink(q1);
  279.     if (qisint(q1) && qisint(q2)) {
  280.         r = qalloc();
  281.         zsub(q1->num, q2->num, &r->num);
  282.         return r;
  283.     }
  284.     q2 = qneg(q2);
  285.     if (qiszero(q1))
  286.         return q2;
  287.     r = qadd(q1, q2);
  288.     qfree(q2);
  289.     return r;
  290. }
  291.  
  292.  
  293. /*
  294.  * Increment a number by one.
  295.  */
  296. NUMBER *
  297. qinc(q)
  298.     NUMBER *q;
  299. {
  300.     NUMBER *r;
  301.  
  302.     r = qalloc();
  303.     if (qisint(q)) {
  304.         zadd(q->num, _one_, &r->num);
  305.         return r;
  306.     }
  307.     zadd(q->num, q->den, &r->num);
  308.     zcopy(q->den, &r->den);
  309.     return r;
  310. }
  311.  
  312.  
  313. /*
  314.  * Decrement a number by one.
  315.  */
  316. NUMBER *
  317. qdec(q)
  318.     NUMBER *q;
  319. {
  320.     NUMBER *r;
  321.  
  322.     r = qalloc();
  323.     if (qisint(q)) {
  324.         zsub(q->num, _one_, &r->num);
  325.         return r;
  326.     }
  327.     zsub(q->num, q->den, &r->num);
  328.     zcopy(q->den, &r->den);
  329.     return r;
  330. }
  331.  
  332.  
  333. /*
  334.  * Add a normal small integer value to an arbitrary number.
  335.  */
  336. NUMBER *
  337. qaddi(q1, n)
  338.     NUMBER *q1;
  339.     long n;
  340. {
  341.     NUMBER addnum;        /* temporary number */
  342.     HALF addval[2];        /* value of small number */
  343.     BOOL neg;        /* TRUE if number is neg */
  344.  
  345.     if (n == 0)
  346.         return qlink(q1);
  347.     if (n == 1)
  348.         return qinc(q1);
  349.     if (n == -1)
  350.         return qdec(q1);
  351.     if (qiszero(q1))
  352.         return itoq(n);
  353.     addnum.num.sign = 0;
  354.     addnum.num.len = 1;
  355.     addnum.num.v = addval;
  356.     addnum.den = _one_;
  357.     neg = (n < 0);
  358.     if (neg)
  359.         n = -n;
  360.     addval[0] = (HALF) n;
  361.     n = (((FULL) n) >> BASEB);
  362.     if (n) {
  363.         addval[1] = (HALF) n;
  364.         addnum.num.len = 2;
  365.     }
  366.     if (neg)
  367.         return qsub(q1, &addnum);
  368.     else
  369.         return qadd(q1, &addnum);
  370. }
  371.  
  372.  
  373. /*
  374.  * Multiply two numbers.
  375.  *    q3 = qmul(q1, q2);
  376.  */
  377. NUMBER *
  378. qmul(q1, q2)
  379.     register NUMBER *q1, *q2;
  380. {
  381.     NUMBER *r;            /* returned value */
  382.     ZVALUE n1, n2, d1, d2;        /* numerators and denominators */
  383.     ZVALUE tmp;
  384.  
  385.     if (qiszero(q1) || qiszero(q2))
  386.         return qlink(&_qzero_);
  387.     if (qisone(q1))
  388.         return qlink(q2);
  389.     if (qisone(q2))
  390.         return qlink(q1);
  391.     if (qisint(q1) && qisint(q2)) {    /* easy results if integers */
  392.         r = qalloc();
  393.         zmul(q1->num, q2->num, &r->num);
  394.         return r;
  395.     }
  396.     n1 = q1->num;
  397.     n2 = q2->num;
  398.     d1 = q1->den;
  399.     d2 = q2->den;
  400.     if (ziszero(d1) || ziszero(d2))
  401.         math_error("Division by zero");
  402.     if (ziszero(n1) || ziszero(n2))
  403.         return qlink(&_qzero_);
  404.     if (!zisunit(n1) && !zisunit(d2)) {    /* possibly reduce */
  405.         zgcd(n1, d2, &tmp);
  406.         if (!zisunit(tmp)) {
  407.             zquo(q1->num, tmp, &n1);
  408.             zquo(q2->den, tmp, &d2);
  409.         }
  410.         zfree(tmp);
  411.     }
  412.     if (!zisunit(n2) && !zisunit(d1)) {    /* again possibly reduce */
  413.         zgcd(n2, d1, &tmp);
  414.         if (!zisunit(tmp)) {
  415.             zquo(q2->num, tmp, &n2);
  416.             zquo(q1->den, tmp, &d1);
  417.         }
  418.         zfree(tmp);
  419.     }
  420.     r = qalloc();
  421.     zmul(n1, n2, &r->num);
  422.     zmul(d1, d2, &r->den);
  423.     if (q1->num.v != n1.v)
  424.         zfree(n1);
  425.     if (q1->den.v != d1.v)
  426.         zfree(d1);
  427.     if (q2->num.v != n2.v)
  428.         zfree(n2);
  429.     if (q2->den.v != d2.v)
  430.         zfree(d2);
  431.     return r;
  432. }
  433.  
  434.  
  435. /*
  436.  * Multiply a number by a small integer.
  437.  *    q2 = qmuli(q1, n);
  438.  */
  439. NUMBER *
  440. qmuli(q, n)
  441.     NUMBER *q;
  442.     long n;
  443. {
  444.     NUMBER *r;
  445.     long d;            /* gcd of multiplier and denominator */
  446.     int sign;
  447.  
  448.     if ((n == 0) || qiszero(q))
  449.         return qlink(&_qzero_);
  450.     if (n == 1)
  451.         return qlink(q);
  452.     r = qalloc();
  453.     if (qisint(q)) {
  454.         zmuli(q->num, n, &r->num);
  455.         return r;
  456.     }
  457.     sign = 1;
  458.     if (n < 0) {
  459.         n = -n;
  460.         sign = -1;
  461.     }
  462.     d = zmodi(q->den, n);
  463.     d = iigcd(d, n);
  464.     zmuli(q->num, (n * sign) / d, &r->num);
  465.     (void) zdivi(q->den, d, &r->den);
  466.     return r;
  467. }
  468.  
  469.  
  470. /*
  471.  * Divide two numbers (as fractions).
  472.  *    q3 = qdiv(q1, q2);
  473.  */
  474. NUMBER *
  475. qdiv(q1, q2)
  476.     register NUMBER *q1, *q2;
  477. {
  478.     NUMBER temp;
  479.  
  480.     if (qiszero(q2))
  481.         math_error("Division by zero");
  482.     if ((q1 == q2) || !qcmp(q1, q2))
  483.         return qlink(&_qone_);
  484.     if (qisone(q1))
  485.         return qinv(q2);
  486.     temp.num = q2->den;
  487.     temp.den = q2->num;
  488.     temp.num.sign = temp.den.sign;
  489.     temp.den.sign = 0;
  490.     temp.links = 1;
  491.     return qmul(q1, &temp);
  492. }
  493.  
  494.  
  495. /*
  496.  * Divide a number by a small integer.
  497.  *    q2 = qdivi(q1, n);
  498.  */
  499. NUMBER *
  500. qdivi(q, n)
  501.     NUMBER *q;
  502.     long n;
  503. {
  504.     NUMBER *r;
  505.     long d;            /* gcd of divisor and numerator */
  506.     int sign;
  507.  
  508.     if (n == 0)
  509.         math_error("Division by zero");
  510.     if ((n == 1) || qiszero(q))
  511.         return qlink(q);
  512.     sign = 1;
  513.     if (n < 0) {
  514.         n = -n;
  515.         sign = -1;
  516.     }
  517.     r = qalloc();
  518.     d = zmodi(q->num, n);
  519.     d = iigcd(d, n);
  520.     (void) zdivi(q->num, d * sign, &r->num);
  521.     zmuli(q->den, n / d, &r->den);
  522.     return r;
  523. }
  524.  
  525.  
  526. /*
  527.  * Return the quotient when one number is divided by another.
  528.  * This works for fractional values also, and in all cases:
  529.  *    qquo(q1, q2) = int(q1 / q2).
  530.  */
  531. NUMBER *
  532. qquo(q1, q2)
  533.     register NUMBER *q1, *q2;
  534. {
  535.     ZVALUE num, den, res;
  536.     NUMBER *q;
  537.  
  538.     if (zisunit(q1->num))
  539.         num = q2->den;
  540.     else if (zisunit(q2->den))
  541.         num = q1->num;
  542.     else
  543.         zmul(q1->num, q2->den, &num);
  544.     if (zisunit(q1->den))
  545.         den = q2->num;
  546.     else if (zisunit(q2->num))
  547.         den = q1->den;
  548.     else
  549.         zmul(q1->den, q2->num, &den);
  550.     zquo(num, den, &res);
  551.     if ((num.v != q2->den.v) && (num.v != q1->num.v))
  552.         zfree(num);
  553.     if ((den.v != q2->num.v) && (den.v != q1->den.v))
  554.         zfree(den);
  555.     if (ziszero(res)) {
  556.         zfree(res);
  557.         return qlink(&_qzero_);
  558.     }
  559.     res.sign = (q1->num.sign != q2->num.sign);
  560.     if (zisunit(res)) {
  561.         q = (res.sign ? &_qnegone_ : &_qone_);
  562.         zfree(res);
  563.         return qlink(q);
  564.     }
  565.     q = qalloc();
  566.     q->num = res;
  567.     return q;
  568. }
  569.  
  570.  
  571. /*
  572.  * Return the absolute value of a number.
  573.  *    q2 = qabs(q1);
  574.  */
  575. NUMBER *
  576. qabs(q)
  577.     register NUMBER *q;
  578. {
  579.     register NUMBER *r;
  580.  
  581.     if (q->num.sign == 0)
  582.         return qlink(q);
  583.     r = qalloc();
  584.     if (!zisunit(q->num))
  585.         zcopy(q->num, &r->num);
  586.     if (!zisunit(q->den))
  587.         zcopy(q->den, &r->den);
  588.     r->num.sign = 0;
  589.     return r;
  590. }
  591.  
  592.  
  593. /*
  594.  * Negate a number.
  595.  *    q2 = qneg(q1);
  596.  */
  597. NUMBER *
  598. qneg(q)
  599.     register NUMBER *q;
  600. {
  601.     register NUMBER *r;
  602.  
  603.     if (qiszero(q))
  604.         return qlink(&_qzero_);
  605.     r = qalloc();
  606.     if (!zisunit(q->num))
  607.         zcopy(q->num, &r->num);
  608.     if (!zisunit(q->den))
  609.         zcopy(q->den, &r->den);
  610.     r->num.sign = !q->num.sign;
  611.     return r;
  612. }
  613.  
  614.  
  615. /*
  616.  * Return the sign of a number (-1, 0 or 1)
  617.  */
  618. NUMBER *
  619. qsign(q)
  620.     NUMBER *q;
  621. {
  622.     if (qiszero(q))
  623.         return qlink(&_qzero_);
  624.     if (qisneg(q))
  625.         return qlink(&_qnegone_);
  626.     return qlink(&_qone_);
  627. }
  628.  
  629.  
  630. /*
  631.  * Invert a number.
  632.  *    q2 = qinv(q1);
  633.  */
  634. NUMBER *
  635. qinv(q)
  636.     register NUMBER *q;
  637. {
  638.     register NUMBER *r;
  639.  
  640.     if (qisunit(q)) {
  641.         r = (qisneg(q) ? &_qnegone_ : &_qone_);
  642.         return qlink(r);
  643.     }
  644.     if (qiszero(q))
  645.         math_error("Division by zero");
  646.     r = qalloc();
  647.     if (!zisunit(q->num))
  648.         zcopy(q->num, &r->den);
  649.     if (!zisunit(q->den))
  650.         zcopy(q->den, &r->num);
  651.     r->num.sign = q->num.sign;
  652.     r->den.sign = 0;
  653.     return r;
  654. }
  655.  
  656.  
  657. /*
  658.  * Return just the numerator of a number.
  659.  *    q2 = qnum(q1);
  660.  */
  661. NUMBER *
  662. qnum(q)
  663.     register NUMBER *q;
  664. {
  665.     register NUMBER *r;
  666.  
  667.     if (qisint(q))
  668.         return qlink(q);
  669.     if (zisunit(q->num)) {
  670.         r = (qisneg(q) ? &_qnegone_ : &_qone_);
  671.         return qlink(r);
  672.     }
  673.     r = qalloc();
  674.     zcopy(q->num, &r->num);
  675.     return r;
  676. }
  677.  
  678.  
  679. /*
  680.  * Return just the denominator of a number.
  681.  *    q2 = qden(q1);
  682.  */
  683. NUMBER *
  684. qden(q)
  685.     register NUMBER *q;
  686. {
  687.     register NUMBER *r;
  688.  
  689.     if (qisint(q))
  690.         return qlink(&_qone_);
  691.     r = qalloc();
  692.     zcopy(q->den, &r->num);
  693.     return r;
  694. }
  695.  
  696.  
  697. /*
  698.  * Return the fractional part of a number.
  699.  *    q2 = qfrac(q1);
  700.  */
  701. NUMBER *
  702. qfrac(q)
  703.     register NUMBER *q;
  704. {
  705.     register NUMBER *r;
  706.     ZVALUE z;
  707.  
  708.     if (qisint(q))
  709.         return qlink(&_qzero_);
  710.     if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
  711.         (q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
  712.             return qlink(q);
  713.     r = qalloc();
  714.     if (qisneg(q)) {
  715.         zmod(q->num, q->den, &z);
  716.         zsub(q->den, z, &r->num);
  717.         zfree(z);
  718.     } else {
  719.         zmod(q->num, q->den, &r->num);
  720.     }
  721.     zcopy(q->den, &r->den);
  722.     r->num.sign = q->num.sign;
  723.     return r;
  724. }
  725.  
  726.  
  727. /*
  728.  * Return the integral part of a number.
  729.  *    q2 = qint(q1);
  730.  */
  731. NUMBER *
  732. qint(q)
  733.     register NUMBER *q;
  734. {
  735.     register NUMBER *r;
  736.  
  737.     if (qisint(q))
  738.         return qlink(q);
  739.     if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
  740.         (q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
  741.             return qlink(&_qzero_);
  742.     r = qalloc();
  743.     zquo(q->num, q->den, &r->num);
  744.     return r;
  745. }
  746.  
  747.  
  748. /*
  749.  * Compute the square of a number.
  750.  */
  751. NUMBER *
  752. qsquare(q)
  753.     register NUMBER *q;
  754. {
  755.     ZVALUE num, den;
  756.  
  757.     if (qiszero(q))
  758.         return qlink(&_qzero_);
  759.     if (qisunit(q))
  760.         return qlink(&_qone_);
  761.     num = q->num;
  762.     den = q->den;
  763.     q = qalloc();
  764.     if (!zisunit(num))
  765.         zsquare(num, &q->num);
  766.     if (!zisunit(den))
  767.         zsquare(den, &q->den);
  768.     return q;
  769. }
  770.  
  771.  
  772. /*
  773.  * Shift an integer by a given number of bits. This multiplies the number
  774.  * by the appropriate power of two.  Positive numbers shift left, negative
  775.  * ones shift right.  Low bits are truncated when shifting right.
  776.  */
  777. NUMBER *
  778. qshift(q, n)
  779.     NUMBER *q;
  780.     long n;
  781. {
  782.     register NUMBER *r;
  783.  
  784.     if (qisfrac(q))
  785.         math_error("Shift of non-integer");
  786.     if (qiszero(q) || (n == 0))
  787.         return qlink(q);
  788.     if (n <= -(q->num.len * BASEB))
  789.         return qlink(&_qzero_);
  790.     r = qalloc();
  791.     zshift(q->num, n, &r->num);
  792.     return r;
  793. }
  794.  
  795.  
  796. /*
  797.  * Scale a number by a power of two, as in:
  798.  *    ans = q * 2^n.
  799.  * This is similar to shifting, except that fractions work.
  800.  */
  801. NUMBER *
  802. qscale(q, pow)
  803.     NUMBER *q;
  804.     long pow;
  805. {
  806.     long numshift, denshift, tmp;
  807.     NUMBER *r;
  808.  
  809.     if (qiszero(q) || (pow == 0))
  810.         return qlink(q);
  811.     if ((pow > 1000000L) || (pow < -1000000L))
  812.         math_error("Very large scale value");
  813.     numshift = zisodd(q->num) ? 0 : zlowbit(q->num);
  814.     denshift = zisodd(q->den) ? 0 : zlowbit(q->den);
  815.     if (pow > 0) {
  816.         tmp = pow;
  817.         if (tmp > denshift)
  818.         tmp = denshift;
  819.         denshift = -tmp;
  820.         numshift = (pow - tmp);
  821.     } else {
  822.         pow = -pow;
  823.         tmp = pow;
  824.         if (tmp > numshift)
  825.             tmp = numshift;
  826.         numshift = -tmp;
  827.         denshift = (pow - tmp);
  828.     }
  829.     r = qalloc();
  830.     if (numshift)
  831.         zshift(q->num, numshift, &r->num);
  832.     else
  833.         zcopy(q->num, &r->num);
  834.     if (denshift)
  835.         zshift(q->den, denshift, &r->den);
  836.     else
  837.         zcopy(q->den, &r->den);
  838.     return r;
  839. }
  840.  
  841.  
  842. /*
  843.  * Return the minimum of two numbers.
  844.  */
  845. NUMBER *
  846. qmin(q1, q2)
  847.     NUMBER *q1, *q2;
  848. {
  849.     if (q1 == q2)
  850.         return qlink(q1);
  851.     if (qrel(q1, q2) > 0)
  852.         q1 = q2;
  853.     return qlink(q1);
  854. }
  855.  
  856.  
  857. /*
  858.  * Return the maximum of two numbers.
  859.  */
  860. NUMBER *
  861. qmax(q1, q2)
  862.     NUMBER *q1, *q2;
  863. {
  864.     if (q1 == q2)
  865.         return qlink(q1);
  866.     if (qrel(q1, q2) < 0)
  867.         q1 = q2;
  868.     return qlink(q1);
  869. }
  870.  
  871.  
  872. /*
  873.  * Perform the logical OR of two integers.
  874.  */
  875. NUMBER *
  876. qor(q1, q2)
  877.     NUMBER *q1, *q2;
  878. {
  879.     register NUMBER *r;
  880.  
  881.     if (qisfrac(q1) || qisfrac(q2))
  882.         math_error("Non-integers for logical or");
  883.     if ((q1 == q2) || qiszero(q2))
  884.         return qlink(q1);
  885.     if (qiszero(q1))
  886.         return qlink(q2);
  887.     r = qalloc();
  888.     zor(q1->num, q2->num, &r->num);
  889.     return r;
  890. }
  891.  
  892.  
  893. /*
  894.  * Perform the logical AND of two integers.
  895.  */
  896. NUMBER *
  897. qand(q1, q2)
  898.     NUMBER *q1, *q2;
  899. {
  900.     register NUMBER *r;
  901.     ZVALUE res;
  902.  
  903.     if (qisfrac(q1) || qisfrac(q2))
  904.         math_error("Non-integers for logical and");
  905.     if (q1 == q2)
  906.         return qlink(q1);
  907.     if (qiszero(q1) || qiszero(q2))
  908.         return qlink(&_qzero_);
  909.     zand(q1->num, q2->num, &res);
  910.     if (ziszero(res)) {
  911.         zfree(res);
  912.         return qlink(&_qzero_);
  913.     }
  914.     r = qalloc();
  915.     r->num = res;
  916.     return r;
  917. }
  918.  
  919.  
  920. /*
  921.  * Perform the logical XOR of two integers.
  922.  */
  923. NUMBER *
  924. qxor(q1, q2)
  925.     NUMBER *q1, *q2;
  926. {
  927.     register NUMBER *r;
  928.     ZVALUE res;
  929.  
  930.     if (qisfrac(q1) || qisfrac(q2))
  931.         math_error("Non-integers for logical xor");
  932.     if (q1 == q2)
  933.         return qlink(&_qzero_);
  934.     if (qiszero(q1))
  935.         return qlink(q2);
  936.     if (qiszero(q2))
  937.         return qlink(q1);
  938.     zxor(q1->num, q2->num, &res);
  939.     if (ziszero(res)) {
  940.         zfree(res);
  941.         return qlink(&_qzero_);
  942.     }
  943.     r = qalloc();
  944.     r->num = res;
  945.     return r;
  946. }
  947.  
  948.  
  949. #if 0
  950. /*
  951.  * Return the number whose binary representation only has the specified
  952.  * bit set (counting from zero).  This thus produces a given power of two.
  953.  */
  954. NUMBER *
  955. qbitvalue(n)
  956.     long n;
  957. {
  958.     register NUMBER *r;
  959.  
  960.     if (n <= 0)
  961.         return qlink(&_qone_);
  962.     r = qalloc();
  963.     zbitvalue(n, &r->num);
  964.     return r;
  965. }
  966.  
  967.  
  968. /*
  969.  * Test to see if the specified bit of a number is on (counted from zero).
  970.  * Returns TRUE if the bit is set, or FALSE if it is not.
  971.  *    i = qbittest(q, n);
  972.  */
  973. BOOL
  974. qbittest(q, n)
  975.     register NUMBER *q;
  976.     long n;
  977. {
  978.     int x, y;
  979.  
  980.     if ((n < 0) || (n >= (q->num.len * BASEB)))
  981.         return FALSE;
  982.     x = q->num.v[n / BASEB];
  983.     y = (1 << (n % BASEB));
  984.     return ((x & y) != 0);
  985. }
  986. #endif
  987.  
  988.  
  989. /*
  990.  * Return the precision of a number (usually for examining an epsilon value).
  991.  * The precision of a number e less than 1 is the positive
  992.  * integer p for which e = 2^-p * f, where 1 <= f < 2.
  993.  * Numbers greater than or equal to one have a precision of zero.
  994.  * For example, the precision of e is 6 if 1/64 <= e < 1/32.
  995.  */
  996. long
  997. qprecision(q)
  998.     NUMBER *q;
  999. {
  1000.     long r;
  1001.  
  1002.     if (qiszero(q) || qisneg(q))
  1003.         math_error("Non-positive number for precision");
  1004.     r = - qilog2(q);
  1005.     return (r < 0 ? 0 : r);
  1006. }
  1007.  
  1008.  
  1009. #if 0
  1010. /*
  1011.  * Return an integer indicating the sign of a number (-1, 0, or 1).
  1012.  *    i = qtst(q);
  1013.  */
  1014. FLAG
  1015. qtest(q)
  1016.     register NUMBER *q;
  1017. {
  1018.     if (!ztest(q->num))
  1019.         return 0;
  1020.     if (q->num.sign)
  1021.         return -1;
  1022.     return 1;
  1023. }
  1024. #endif
  1025.  
  1026.  
  1027. /*
  1028.  * Determine whether or not one number exactly divides another one.
  1029.  * Returns TRUE if the first number is an integer multiple of the second one.
  1030.  */
  1031. BOOL
  1032. qdivides(q1, q2)
  1033.     NUMBER *q1, *q2;
  1034. {
  1035.     if (qiszero(q1))
  1036.         return TRUE;
  1037.     if (qisint(q1) && qisint(q2)) {
  1038.         if (qisunit(q2))
  1039.             return TRUE;
  1040.         return zdivides(q1->num, q2->num);
  1041.     }
  1042.     return zdivides(q1->num, q2->num) && zdivides(q2->den, q1->den);
  1043. }
  1044.  
  1045.  
  1046. /*
  1047.  * Compare two numbers and return an integer indicating their relative size.
  1048.  *    i = qrel(q1, q2);
  1049.  */
  1050. FLAG
  1051. qrel(q1, q2)
  1052.     register NUMBER *q1, *q2;
  1053. {
  1054.     ZVALUE z1, z2;
  1055.     long wc1, wc2;
  1056.     int sign;
  1057.     int z1f = 0, z2f = 0;
  1058.  
  1059.     if (q1 == q2)
  1060.         return 0;
  1061.     sign = q2->num.sign - q1->num.sign;
  1062.     if (sign)
  1063.         return sign;
  1064.     if (qiszero(q2))
  1065.         return !qiszero(q1);
  1066.     if (qiszero(q1))
  1067.         return -1;
  1068.     /*
  1069.      * Make a quick comparison by calculating the number of words resulting as
  1070.      * if we multiplied through by the denominators, and then comparing the
  1071.      * word counts.
  1072.      */
  1073.     sign = 1;
  1074.     if (qisneg(q1))
  1075.         sign = -1;
  1076.     wc1 = q1->num.len + q2->den.len;
  1077.     wc2 = q2->num.len + q1->den.len;
  1078.     if (wc1 < wc2 - 1)
  1079.         return -sign;
  1080.     if (wc2 < wc1 - 1)
  1081.         return sign;
  1082.     /*
  1083.      * Quick check failed, must actually do the full comparison.
  1084.      */
  1085.     if (zisunit(q2->den))
  1086.         z1 = q1->num;
  1087.     else if (zisone(q1->num))
  1088.         z1 = q2->den;
  1089.     else {
  1090.         z1f = 1;
  1091.         zmul(q1->num, q2->den, &z1);
  1092.     }
  1093.     if (zisunit(q1->den))
  1094.         z2 = q2->num;
  1095.     else if (zisone(q2->num))
  1096.         z2 = q1->den;
  1097.     else {
  1098.         z2f = 1;
  1099.         zmul(q2->num, q1->den, &z2);
  1100.     }
  1101.     sign = zrel(z1, z2);
  1102.     if (z1f)
  1103.         zfree(z1);
  1104.     if (z2f)
  1105.         zfree(z2);
  1106.     return sign;
  1107. }
  1108.  
  1109.  
  1110. /*
  1111.  * Compare two numbers to see if they are equal.
  1112.  * This differs from qrel in that the numbers are not ordered.
  1113.  * Returns TRUE if they differ.
  1114.  */
  1115. BOOL
  1116. qcmp(q1, q2)
  1117.     register NUMBER *q1, *q2;
  1118. {
  1119.     if (q1 == q2)
  1120.         return FALSE;
  1121.     if ((q1->num.sign != q2->num.sign) || (q1->num.len != q2->num.len) ||
  1122.         (q2->den.len != q2->den.len) || (*q1->num.v != *q2->num.v) ||
  1123.         (*q1->den.v != *q2->den.v))
  1124.             return TRUE;
  1125.     if (zcmp(q1->num, q2->num))
  1126.         return TRUE;
  1127.     if (qisint(q1))
  1128.         return FALSE;
  1129.     return zcmp(q1->den, q2->den);
  1130. }
  1131.  
  1132.  
  1133. /*
  1134.  * Compare a number against a normal small integer.
  1135.  * Returns 1, 0, or -1, according to whether the first number is greater,
  1136.  * equal, or less than the second number.
  1137.  *    n = qreli(q, n);
  1138.  */
  1139. FLAG
  1140. qreli(q, n)
  1141.     NUMBER *q;
  1142.     long n;
  1143. {
  1144.     int sign;
  1145.     ZVALUE num;
  1146.     HALF h2[2];
  1147.     NUMBER q2;
  1148.  
  1149.     sign = ztest(q->num);        /* do trivial sign checks */
  1150.     if (sign == 0) {
  1151.         if (n > 0)
  1152.             return -1;
  1153.         return (n < 0);
  1154.     }
  1155.     if ((sign < 0) && (n >= 0))
  1156.         return -1;
  1157.     if ((sign > 0) && (n <= 0))
  1158.         return 1;
  1159.     n *= sign;
  1160.     if (n == 1) {            /* quick check against 1 or -1 */
  1161.         num = q->num;
  1162.         num.sign = 0;
  1163.         return (sign * zrel(num, q->den));
  1164.     }
  1165.     num.sign = (sign < 0);
  1166.     num.len = 1 + (n >= BASE);
  1167.     num.v = h2;
  1168.     h2[0] = (n & BASE1);
  1169.     h2[1] = (n >> BASEB);
  1170.     if (zisunit(q->den))    /* integer compare if no denominator */
  1171.         return zrel(q->num, num);
  1172.     q2.num = num;
  1173.     q2.den = _one_;
  1174.     q2.links = 1;
  1175.     return qrel(q, &q2);    /* full fractional compare */
  1176. }
  1177.  
  1178.  
  1179. /*
  1180.  * Compare a number against a small integer to see if they are equal.
  1181.  * Returns TRUE if they differ.
  1182.  */
  1183. BOOL
  1184. qcmpi(q, n)
  1185.     NUMBER *q;
  1186.     long n;
  1187. {
  1188.     long len;
  1189.  
  1190.     len = q->num.len;
  1191.     if ((len > 2) || qisfrac(q) || (q->num.sign != (n < 0)))
  1192.         return TRUE;
  1193.     if (n < 0)
  1194.         n = -n;
  1195.     if (((HALF)(n)) != q->num.v[0])
  1196.         return TRUE;
  1197.     n = ((FULL) n) >> BASEB;
  1198.     return (((n != 0) != (len == 2)) || (n != q->num.v[1]));
  1199. }
  1200.  
  1201.  
  1202. /*
  1203.  * Number node allocation routines
  1204.  */
  1205.  
  1206. #define    NNALLOC    1000
  1207.  
  1208. union allocNode {
  1209.     NUMBER    num;
  1210.     union allocNode    *link;
  1211. };
  1212.  
  1213. static union allocNode    *freeNum;
  1214.  
  1215.  
  1216. NUMBER *
  1217. qalloc()
  1218. {
  1219.     register union allocNode *temp;
  1220.  
  1221.     if (freeNum == NULL) {
  1222.         freeNum = (union allocNode *)
  1223.         malloc(sizeof (NUMBER) * NNALLOC);
  1224.         if (freeNum == NULL)
  1225.             math_error("Not enough memory");
  1226.         freeNum[NNALLOC-1].link = NULL;
  1227.         for (temp=freeNum+NNALLOC-2; temp >= freeNum; --temp) {
  1228.             temp->link = temp+1;
  1229.         }
  1230.     }
  1231.     temp = freeNum;
  1232.     freeNum = temp->link;
  1233.     temp->num.links = 1;
  1234.     temp->num.num = _one_;
  1235.     temp->num.den = _one_;
  1236.     return &temp->num;
  1237. }
  1238.  
  1239.  
  1240. void
  1241. qfreenum(q)
  1242.     register NUMBER *q;
  1243. {
  1244.     union allocNode *a;
  1245.  
  1246.     if (q == NULL)
  1247.         return;
  1248.     zfree(q->num);
  1249.     zfree(q->den);
  1250.     a = (union allocNode *) q;
  1251.     a->link = freeNum;
  1252.     freeNum = a;
  1253. }
  1254.  
  1255. /* END CODE */
  1256.